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Abstract 

The quantum dynamics of a one-dimensional bosonic Josephson junction is studied by solv- 
ing the time-dependent many-boson Schrodinger equation numerically exactly. Already for weak 
interparticle interactions and on short time scales, the commonly-employed mean- field and many- 
body methods are found to deviate substantially from the exact dynamics. The system exhibits 
rich many-body dynamics like enhanced tunneling and a novel equilibration phenomenon of the 
junction depending on the interaction, attributed to a quick loss of coherence. 
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Recent experiments on interacting Bose-Einstein condensates in double-well traps have 
led to some of the most exciting results in quantum physics, including matter-wave inter- 
'erometry ll, |2|], squeezing and entanglement 31,1^ as well as work on high- precision sensors 
sl- Particular attention has been paid to tunneling phenomena of interacting Bose-Einstein 
condensates in double-wells, which in this context are referred to as bosonic Josephson junc- 



tions. Explicitly, Josephson oscillations and sel 
3ose- Einstein condensates have been predicted |6 
s], Q], drawing intensive interest, see, e.g., 10, 11 



-trapping (suppression of tunneling) with 
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and recently realized in experiments 
1^ . [isl and references therein. 



For the first time in literature we provide the numerically exact many-body quantum 
dynamics of a one-dimensional (ID) bosonic Josephson junction in this work. This is made 
possible by a breakthrough in the solution of the time- dependent many-boson Schrodinger 
equation. We use the exact solution to check the current understanding of bosonic Josephson 
junctions - commonly described by the popular Gross-Pitaevskii (GP) mean-field theory 
and the Bose-Hubbard (BH) many-body model - and to find novel phenomena. The results 
of the GP and BH theories are found to deviate substantially from the full many-body 
solution, already for weak interactions and on short time scales. In particular, the well- 
known self-trapping effect is greatly reduced. We attribute these findings to a quick loss of 
the junction's coherence not captured by the common methods. For stronger interactions 
and on longer time scales, we find a novel equilibration dynamics in which the density and 
other observables of the junction tend towards stationary values. We show that the dynamics 
of bosonic Josephson junctions is much richer than what is currently known. 

To compute the time evolution of the system, we solve the time-dependent many-boson 
Schrodinger equation by using the multiconfigurational time-dependent Hartree for bosons 
(MCTDHB) method [id]. In the MCTDHB(M) method the time-dependent many-boson 
wavefunction is expanded in all time- dependent permanents generated by distribut- 

ing bosons over M time-dependent orbitals {0j(x, t)}. n = {ni,n2, ■ ■ ■ um) collects the 
occupation numbers. The MCTDHB wavefunction thus reads {"^{t)) = Xln '^^^ 
expansion coefficients {Cft{t)} and the orbitals {4>i{x, t)} are determined by the Dirac-Frenkel 
time-dependent variational principle |l^. The present results are obtained by using a novel 
mapping of the many-boson configuration space in combination with a parallel implementa- 
tion of MCTDHB, allowing the efficient handling of millions of time-dependent, optimized 
permanents We note that GP theory is contained in the MCTDHB framework as the 
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special case M = 1. The many-body wavefunction then becomes a single permanent and the 
dynamics is restricted to remaining condensed at all times. The BH model for this system 
employs two orbitals and is thus, in principle, capable of describing correlations. However, 
as we shall see below, the fact that the BH orbitals are time-independent causes the BH 
model to underestimate correlations. 

With the time-dependent many-boson wavefunction at hand, any quantity of inter- 

est of the interacting many-boson system can be computed. Here we focus on the evolution 
of the following quantities to analyze the dynamics of the Josephson junction. The reduced 
one-body density matrix of the system is defined by p^^\x\x';t) = (\l/(t)|^'^(a;')^(x)|\E'(i(:)), 
where ^(x) is the usual bosonic field operator annihilating a particle at position x. Its 
diagonal part, p{x,t) = p^-^\x\x' = x;t), is simply the density of the system. As is 
common in the analysis of bosonic Josephson junctions, the "survival probability" of the 
system in, e.g., the left well, is obtained by integrating the density over the left well, 
piit) = p{x,t)dx. Furthermore, the eigenvalues n[^^ of p^^\x\x';t) determine the 

extent to which the system is condensed or fragmented 18|, ll9l] . Finally, the first-order cor- 



^^Hx', x; t) = p'--^\x\x'; t)/ a/ p(x, t)p{x', t) quantifies the system's degree of 



spatial coherence 



9«( 


x', X 


20, 


21|. 



We now turn to the details of the ID bosonic Josephson junction considered in this work. 
It is convenient to use dimensionless units defined by dividing the Hamiltonian by 
where m is the mass of a boson, e.g., ^^Rb and L is a length scale, e.g., L = 1pm. One 
unit of energy then corresponds to 116 Hz. We have a ID realization of an experimental 
setup similar to the one in Ref. 8|] in mind. The full many-body Hamiltonian then reads 
H = X^ili h{xi) + X]i<i ^i^i ~ ^j)-> where h{x) = + V{x), with a trapping potential 

V{x) and an interparticle interaction potential W{x — x') = Xq6{x — x'). 

The double-well potential V{x) is generated by connecting two harmonic potentials 
V±(x) = ± 2)^ with a cubic spline in the region < 0.5. This results in a symmet- 
ric double-well potential with barrier height V{0) = 1.667. The lowest four single-particle 
energy levels ei = 0.473,62 = 0.518,63 = 1.352 and 64 = 1.611 of V{x) are lower than the 
barrier. 

Left- and right-localized orbitals (f)L,R can be constructed from the single-particle ground 
state and the first excited state of V{x). (pL and (pR determine the parameters U = Xq f \4>l\'^, 
J = — f (j)*j;^h(f)R, the Rabi oscillation period t^abi = vr/J and the often employed interaction 



3 



parameters, A = UN/ (2 J) and U / J . In this work we use the interaction parameter 
A = Ao(A^ — 1), which appears naturally in the full many-body treatment, and quote the 
corresponding values for A and U / J . Within the framework of two-mode GP theory, a state, 
which is initially localized in one well, is predicted to remain s elf -trapped ii A > A^, = 2 6|, |7|. 
We will consider interaction strengths below, in the vicinity of and above Ac. 

In all our computations the system is prepared at t = as the many-body ground state of 
the potential V+{x) and then propagated in the potential V{x). Within the BH framework 
this procedure amounts to starting from the state in which all bosons occupy the orbital 0^. 

We begin our studies with a weak interaction strength A = 0.152, leading to U/J = 0.140 
(0.027) and A = 1.40 (1.35) for = 20 (100) bosons, which is well below the transition 
point for self-trapping A^ = 2. In the upper two panels of Fig. [T] the full many-body 
(solid blue lines) results for pL{t) are shown together with those of GP (solid black lines) 
and BH (solid magenta lines) theory. The full many-body dynamics is governed by three 
different time scales. On a time scale of the order of a Rabi cycle, piit) performs large- 
amplitude oscillations about pl = 0.5, the long time average of piif)- The amphtude of 
these oscillations is damped out on a time scale of a few Rabi cycles and marks the beginning 
of a collapse and revival (not shown) sequence j^, confirmed here on the full many-body 
level. On top of these slow large-amplitude oscillations, a higher frequency with a small 
amplitude can be seen. In a single-particle picture these high frequency oscillations can be 
related to contributions from higher excited states in the initial wavefunction. However, 
a single particle picture fails to describe the dynamics, as we shall now show. While the 
initial wavefunction |\I'(t = 0)) is practically condensed - the fragmentation of the system 
is less than 10~^ (10~^) for = 20 (100) bosons - the propagated wavefunction |^E'(t)) 
quickly becomes fragmented. The fragmentation increases to about 33% (26%) at t = Stf^abi 
for = 20 (100) particles, making a many-body treatment indispensable, already at this 
weak interaction strength. The respective GP results (solid black lines) oscillate back and 
forth at a frequency close to the Rabi frequency and resemble the full many-body dynamics 
only on a time scale shorter than half a Rabi cycle. The poor quality of the GP mean-field 
approximation is, of course, due to the fact that the exact wavefunction starts to fragment 
while the GP dynamics remains condensed by construction. 

The BH (solid magenta line) result for piit) reproduces many features of that of the full 
many-body solution at this interaction strength for both A^ = 20 and A^ = 100 particles. 
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The large-amplitude oscillations collapse over a period of a few Rabi cycles and revive at a 
later stage (not shown). Also the BH solution quickly becomes fragmented, starting from the 
left localized state, which is totally condensed. The fragmentation of the BH wavefunction 
for = 20 (100) particles at t = St^au is essentially the same as the respective value of the 
exact solution. However, differences between the exact and the BH result are visible even 
on time scales less than half a Rabi cycle. Not only are the amplitudes obviously different, 
but also the frequencies contained in pL{t)- Furthermore, the BH solutions do not exhibit a 
high frequency oscillation on top of the slow large-amplitude oscillations; a difference which 
is related to the fact that the BH orbitals are time-independent and thus, not determined 
variationally at each point in time. Note that Piit) is a quantity in which all spatial degrees 
of freedom have been integrated out. Visible differences in PL{t) imply that it is not only 
the densities p{x, t) which must differ, but also all correlation functions. 

The insets of Fig. [T](a),(b) demonstrate the convergence of the many-body dynamics 
results. In particular and somewhat unexpectedly, the number of time- dependent orbitals 
needed to describe the bosonic Josephson junction dynamics quantitatively is M = 4, even 
below the transition point for self-trapping. These orbitals are determined variationally at 
each point in time, implying that any method using time-independent orbitals will need 
substantially more orbitals to achieve the same accuracy. 

One of the central phenomena often discussed in the context of bosonic Josephson junc- 
tions is the celebrated transition to self-trapping . In what follows we would like 
to study the dynamics of a bosonic Josephson junction in the self-trapping regime from the 
full many-body perspective. 

The interaction strength is taken to be A = 0.245, leading to U/J = 0.226 (0.043) and 
A = 2.26 (2.17) for A^ = 20 (100)_particles. Hence, the system is just above the critical 
value for self-trapping Ac = 2 0, IJ. The results for A^ = 20 and A^ = 100 are collected 
in Fig. [l](c),(d). We find that the full many-body solutions (solid blue lines) exhibit indeed 
some self-trapping on the time scale shown. The fragmentation of the condensate for A^ = 20 
(100) bosons increases from initially less than 10"'^ (10^^) to about 28% (18%) after three 
Rabi cycles. Note that the system is now less fragmented than for weaker interactions 
after the same period of time. Nevertheless, GP (solid black lines) theory is - as before - 
inapplicable, even on time scales shorter than tjiabif^- The BH (solid magenta lines) results 
deviate from the true dynamics even earlier. They greatly overestimate the self-trapping 
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and coherence of the condensate. According to the BH model the condensate would only be 
13% (11%) fragmented for = 20 (100) at t = StRabi, which is not the case. This trend also 
continues for stronger interactions, see below. The following general statement about the 
relationship between self-trapping and coherence can be inferred from our full many-body 
results: Self-trapping is only present as long as the system remains coherent. We find this 
statement to be true at all interaction strengths and all particle numbers considered in this 
work. 

We now turn to the case of stronger interactions, A = 4.9, which is well above the 
self-trapping transition point. This leads to U/J = 9.55 (0.869) and A = 47.8 (43.4) for 
A^ = 10 (100) bosons. Note that we now use ten instead of twenty bosons to demonstrate 
convergence. The energy per particle of the full many-body wavefunction is now E/N = 1.22 
(1.28) for A^ = 10 (100) bosons, which is still below the barrier height ^(0) = 1.667. 
What do we expect to happen? According to two-mode GP theory, the density should 
remain trapped in the initial well for any interaction strength A ^ Ac 6|, LZl]. Similarly, the 
BH model predicts ever increasing tunneling times since the left- ancl_right-localized states 
become eigenstates of the BH Hamiltonian in the limit f//J — oo [6]. These predictions 
are incorrect for stronger interactions since the repulsive interaction leads to fragmentation 
and broadening of the initial wavefunction, thereby facilitating tunneling. The impact of 
the potential barrier therefore decreases with increasing interaction strength. We will now 
show that a very intricate dynamics results. 

Fig. [2](top) shows the full many-body results for A^ = 10 (solid blue line) and A^ = 100 
(solid green line) bosons together with those of the BH (solid magenta line) model. The two 
BH results lie on top of each other. In complete contrast to the BH dynamics, for which 
Piit) remains trapped in the left well, the full many-body dynamics shows no self-trapping. 
Instead, an equilibration phenomenon emerges, in which the density of the system tends to 
be equally distributed over both wells. 

The system's full many-body dynamics is again strongly fragmented as can be seen in 
Fig. [2](bottom), which depicts the natural-orbital occupations nf^ (solid blue lines) for 
A^ = 10 particles. This rules out any description of the system by GP mean-field theory. 
Also shown are the natural-orbital occupations of the BH (solid magenta lines) model, which 
wrongly describes a fully condensed system. 

The strong fragmentation of the system implies the presence of strong correlations. This 
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can be seen in the two upper panels of Fig. [3l which show the full many-body result for the 
first-order correlation function g^^\x',X]t) of = 10 bosons at times t = (top left) and 
t = lOtRabi (top right). The fragmentation of the initial state is only ^ 2%, leading to an 
almost fiat g^^\x',x;0). This reflects the fact that the system is initially coherent over its 
entire extent. At t = lOtjiabi the coherence of the system is completely lost even on length 
scales much shorter than its size, see upper right panel of Fig. [31 Note that also g^^\x', x; t) 
tends to equilibrate. The respective BH results for g^-^\x',x;t) are shown in the two lower 
panels of Fig. [3] and in contrast display no visible loss of coherence. 

Let us briefly summarize. We have obtained exact results for the full many-body dynamics 
of ID bosonic Josephson junctions. The dynamics is found to be much richer than previously 
reported. In particular, the predictions of the commonly-employed Gross-Pitaevskii and 
Bose-Hubbard theories are found to differ substantially from the exact results, already after 
short times and relatively weak interactions. These differences are associated with the 
development of fragmentation and correlations not captured by the standard theories. For 
stronger interactions, where the standard theories predict coherence and self-trapping, we 
find a completely different dynamics. The system becomes fragmented, spatial coherence is 
lost and a long-time equilibration of the junction emerges. We hope our results stimulate 
experiments. 
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FIG. 1: (Color online) Full quantum dynamics of a ID bosonic Josephson junction below and 
above the transition to self-trapping. Shown is the full many-body result (solid blue lines) for the 
probability of finding a boson in the left well, PL{t)- For comparison, the respective GP (solid 
black lines) and BH (solid magenta lines) results are shown as well. The parameter values are: (a) 
AT = 20, A = 0.152 and (b) N = 100, A = 0.152 (below the self-trapping transition), (c) N = 20, 
A = 0.245 and (d) N = 100, A = 0.245 (above the self-trapping transition). The GP and BH 
results are found to deviate from the full many-body results already after short times. The insets 
show the convergence of the full many-body results. (a),(c): M = 2 (solid purple line), M = 4 
(solid red line), M = 6 (solid green line), M = 8 (solid blue line). The M = 2 results are seen to 
deviate slightly from the converged results for M > 4. (b),(d): The results for M = 2 (solid purple 
line) and M = 4 (solid blue line) are shown. All quantities shown are dimensionless. 
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FIG. 2: (Color online) Emergence of equilibration of the density at interaction strength A = 4.9. 
Top: same as Fig.[Tl but for = 10 (solid blue line) and N = 100 (solid green line). The respective 
BH (solid magenta lines) results are on top of each other. In contrast to the BH dynamics which 
is completely self-trapped, the full many-body dynamics is not. PL{t) tends towards its long- 
time average pi = 0.5. For N = 100 particles M = 4 orbitals were used. The inset shows the 
convergence of the full many-body solution for = 10 bosons: M = 4 (solid black line), M = 10 
(solid blue line), M = 12 (solid red line). The M = 4 result follows the trend of the converged 
M = 12 result. Bottom: corresponding natural orbital occupations for = 10 bosons. The system 
becomes fragmented and roughly four natural orbitals are macroscopically occupied. All quantities 
are dimensionless. 
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FIG. 3: (Color online) Dynamics of the first order correlation function for A = 4.9 at which the 
equilibration phenomenon of Fig. [2] occurs. Shown is \g'^'^\x' , x; t)p of = 10 bosons at different 
times. Top left: full many-body result at t = 0. The initial state exhibits coherence over the entire 
extent of the system. Top right: full many-body result at t = lOtfiau- The coherence is lost. The 
system is incoherent even on short length scales. Bottom left: BH result at t = 0. Bottom right: 
BH result at t = Wtjiabi- In contrast to the full many-body result, the BH wavefunction remains 
completely coherent. 
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